Disclaimer

This work is very preliminary as I get back into the coding swing of things. Data wrangling and figure generation will be done via R, but the rest of the project will be done using good ol’ microsoft products. This is just an entry point into data crunching and should by no means be considered a final product.

Also, I’m not great at this but whatever. I could automate this, but I’ll figure that out shortly!

Need to update figures & analyze sites for seasonal trends

Methodology

SNOTEL data was provided by the NRCS. Data was cleaned by removing outliers that are likely implausible; any year with more than 15 observations missing was removed. Temperatures were adjusted using the Morrisey method for stations identified by Ma et al (2019) due to SNOTEL temperature sensor changes, with the adjustment applied to pre-sensor change data. Daily mean observations were detrended to determine whether values were increasing or decreasing from the entire time series trend. Daily mean temperatures were first averaged by water year, with all water year means then averaged by day of water year. The mean temperature by day for the period of record was averaged. To find the standard deviation, the daily mean temperatures by water year was subtracted from the averaged mean temperature by day for the period of record. All water year means averaged by day of water year were subtracted from the temperature mean. The resulting values were then added together to find the “residual” of the daily mean temperatures by water year. The standard deviation was then computed from those residuals, with trends analyzed by Mann‐Kendall significance test and Theil‐Sen’s rate of change. Significant trends are identified with p-values of less than 0.10.

Morrisey Method

The Morrisey Method is taken from Ma, Fassnacht and Kampf..

In R script: T(adjusted) = 5.3x10(-7)xT(old)4+3.72x10(-5)xT(old)3-2.16x10(-3)xT(old)2-7.32x10^(-2)xT(old)+1.37

In the Ma et al. spreadsheet, H1 is Morrisey, H2 is Oiler

San Juan Area SNOTEL sites:

Beartown 327 Original 4/18/2005

Cascade #2 389 Morrisey 6/18/2004

Cumbres Trestle 431 Oyler -> Morrisey 6/8/2005

Idarado 538 Morrisey 7/12/2004

Lone Cone 589 Oyler -> Morrisey 6/22/2005

Middle Creek 624 Morrisey 6/26/2006

Mineral Creek 629 Oyler ?? -> Morrisey 6/22/2004

Molas Lake 632 NSCE-Morrisey, Bias-Original 10/2/2003

Red Mountain Pass 713 Morrisey 8/18/2004

Scotch Creek 739 Morrisey 6/15/2004

Slumgullion 762 Morrisey 6/26/2006

Spud Mountain 780 Morrisey 6/28/2004

Stump Lakes 797 Morrisey 7/22/2005

Upper Rio Grande 839 Oyler -> Morrisey 5/26/2004 *Why is this in red?

Upper San Juan 840 Morrisey 4/7/2004

Vallecito 843 Morrisey 7/22/2005

Wolf Creek Summit 874 Morrisey 7/12/2004

Data from thesis research:


SNOTEL_san_juan_Area <- snotel_download(site_id = c(327, 389, 431, 538, 589, 624, 629, 632, 713, 739, 762, 780, 797, 839, 840, 843, 874), path = tempdir('../data'), internal = TRUE)

write.csv(SNOTEL_san_juan_Area,"C:/Users/13074/Documents/ESS580/thesis_project/San_Juan_area/data_raw/snotel_san_juans.csv", row.names = FALSE) #write in the raw data

Beartown 327

Original

SNOTEL_san_juan_Area <- read.csv("C:/Users/13074/Documents/ESS580/thesis_project/San_Juan_area/data_raw/snotel_san_juans.csv", header = TRUE)

snotel_327 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "327")
#str(snotel_327) # check the date, usually a character.  

snotel_327$Date <- as.Date(snotel_327$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_327_clean <- snotel_327 %>% # filter for the timeframe
  filter(Date >= "1983-08-10" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_327_clean <- snotel_327_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_327_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_327_clean <- snotel_327_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_min > -40)
ggplot(snotel_327_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
snotel_327_cull_count <- snotel_327_clean %>% 
  filter(temperature_min < -40) %>% 
  count(waterYear)

snotel_327_cull_count
## # A tibble: 0 x 2
## # Groups:   waterYear [0]
## # ... with 2 variables: waterYear <dbl>, n <int>
# filtering for too few observations in a year
snotel_327_cull_count_days <- snotel_327_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_327_cull_count_days
## # A tibble: 3 x 2
## # Groups:   waterYear [3]
##   waterYear     n
##       <dbl> <int>
## 1      1983    52
## 2      1994   337
## 3      2003   346

1983, 1994, 2003 need to be culled.

snotel_327_clean_culled <- snotel_327_clean %>% 
  filter(waterYear != "1983" & waterYear != "1994" & waterYear != "2003") #%>% 
  #filter(temperature_mean > -49)

Beartown (327) NOT ADJUSTED.

2005-04-18 installed ext range sensor. Original

snotel_327_adjusted <- snotel_327_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2005-04-18", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

327 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_327 <- snotel_327_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_327 <- yearly_wy_aver_327 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_327 <- daily_wy_aver_327 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_327$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_327 <-daily_wy_aver_327 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_327$date_temp <- signif(daily_wy_aver2_327$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_327, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

327 SD

standard_dev_327 <- daily_wy_aver_327 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_327 <- standard_dev_327 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_327 <- standard_dev_all_327 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_327 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 8.215357
1985 8.122497
1986 7.129599
1987 7.619777
1988 8.257318
1989 8.044976
1990 7.675963
1991 7.684970
1992 7.079517
1993 7.642300
1995 7.304999
1996 7.681058
1997 7.711950
1998 8.105700
1999 6.825322
2000 7.526413
2001 7.912097
2002 8.099855
2004 7.622635
2005 7.225545
2006 7.107924
2007 7.426444
2008 7.802192
2009 6.893765
2010 7.820347
2011 7.632307
2012 7.466342
2013 8.000100
2014 7.260505
2015 6.671562
2016 7.474460
2017 7.072348
2018 6.933787
2019 7.658054
2020 7.769482
2021 7.841403
2022 7.318492
ggplot(standard_dev_all_327, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 327 average temperatures for water years 2005-2021

MK & SS for 327 (non-corrected)

sd_mk_327 <- mk.test(standard_dev_all_327$sd_2)
print(sd_mk_327)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_327$sd_2
## z = -1.818, n = 37, p-value = 0.06907
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -140.0000000 5846.0000000   -0.2102102
sd_sens_327 <- sens.slope(standard_dev_all_327$sd_2)
print(sd_sens_327)
## 
##  Sen's slope
## 
## data:  standard_dev_all_327$sd_2
## z = -1.818, n = 37, p-value = 0.06907
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.025461316  0.001275883
## sample estimates:
## Sen's slope 
## -0.01223021

Cascade #2 389

Morrisey 6/18/2004

snotel_389 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "389")
#str(snotel_389) # check the date, usually a character.  

snotel_389$Date <- as.Date(snotel_389$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_389_clean <- snotel_389 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_389_clean <- snotel_389_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_389_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_389_clean <- snotel_389_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_min > -40) %>% 
  filter(temperature_mean > -40)
ggplot(snotel_389_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

temp_389_xts <- xts(snotel_389_clean$temperature_mean, order.by = snotel_389_clean$Date)

dygraph(temp_389_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
snotel_389_cull_count <- snotel_389_clean %>% 
  filter(temperature_min > -40) %>% 
  count(waterYear)

snotel_389_cull_count
## # A tibble: 40 x 2
## # Groups:   waterYear [40]
##    waterYear     n
##        <dbl> <int>
##  1      1983   325
##  2      1984   350
##  3      1985   350
##  4      1986   365
##  5      1987   358
##  6      1988   366
##  7      1989   330
##  8      1990   348
##  9      1991   364
## 10      1992   366
## # ... with 30 more rows
# filtering for too few observations in a year
snotel_389_cull_count_days <- snotel_389_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_389_cull_count_days
## # A tibble: 9 x 2
## # Groups:   waterYear [9]
##   waterYear     n
##       <dbl> <int>
## 1      1983   325
## 2      1989   330
## 3      1990   348
## 4      1994   310
## 5      2006   330
## 6      2010   334
## 7      2014   314
## 8      2015   185
## 9      2016   334

1983, 1989, 1990, 1994, 2006, 2010, 2014, 2015, 2016 need to be culled.

snotel_389_clean_culled <- snotel_389_clean %>% 
  filter(waterYear != "1983" & waterYear != "1990" & waterYear != "1994" & waterYear != "2006" & waterYear != "2010" & waterYear != "2014" & waterYear != "2015" & waterYear != "2016") #%>% 
  #filter(temperature_mean > -49)

Cascade #2 (389) NOT ADJUSTED.

Morrisey 6/18/2004

snotel_389_adjusted <- snotel_389_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2004-06-18", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

389 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_389 <- snotel_389_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_389 <- yearly_wy_aver_389 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_389 <- daily_wy_aver_389 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_389$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_389 <-daily_wy_aver_389 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_389$date_temp <- signif(daily_wy_aver2_389$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_389, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

389 SD

standard_dev_389 <- daily_wy_aver_389 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_389 <- standard_dev_389 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_389 <- standard_dev_all_389 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_389 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 10.776321
1985 7.879611
1986 10.887499
1987 9.408888
1988 10.711609
1989 10.418326
1991 10.101153
1992 8.637142
1993 8.884440
1995 9.371322
1996 9.429341
1997 10.027003
1998 9.698491
1999 9.029327
2000 9.288091
2001 10.630687
2002 10.346726
2003 9.727202
2004 9.282545
2005 8.472227
2007 9.686674
2008 9.541070
2009 8.659960
2011 9.334867
2012 9.418612
2013 9.751020
2017 9.144382
2018 8.911243
2019 9.427120
2020 9.880776
2021 9.863539
2022 9.441803
ggplot(standard_dev_all_389, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 389 average temperatures for water years 2005-2021

MK & SS for 389 (non-corrected)

sd_mk_389 <- mk.test(standard_dev_all_389$sd_2)
print(sd_mk_389)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_389$sd_2
## z = -0.95677, n = 32, p-value = 0.3387
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  -60.0000000 3802.6666667   -0.1209677
sd_sens_389 <- sens.slope(standard_dev_all_389$sd_2)
print(sd_sens_389)
## 
##  Sen's slope
## 
## data:  standard_dev_all_389$sd_2
## z = -0.95677, n = 32, p-value = 0.3387
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.04702986  0.01336703
## sample estimates:
## Sen's slope 
## -0.01646233

Corrected

389 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_389_ad <- snotel_389_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_389_ad <- yearly_wy_aver_389_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_389_ad <- daily_wy_aver_389_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_389_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_389_ad <-daily_wy_aver_389_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_389_ad$date_temp_ad <- signif(daily_wy_aver2_389_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_389_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

389 SS (corrected)

standard_dev_389_ad <- daily_wy_aver_389_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_389_ad <- standard_dev_389_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_389_ad <- standard_dev_all_389_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_389_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 10.061639
1985 7.175032
1986 9.930924
1987 8.682951
1988 9.868631
1989 9.733860
1991 9.472969
1992 8.002025
1993 8.283300
1995 8.657379
1996 8.697682
1997 9.311693
1998 8.929771
1999 8.328247
2000 8.510889
2001 9.808164
2002 9.563945
2003 8.940614
2004 8.587303
2005 8.472294
2007 9.686174
2008 9.541322
2009 8.659292
2011 9.334237
2012 9.418205
2013 9.750812
2017 9.144257
2018 8.911176
2019 9.426664
2020 9.880403
2021 9.863015
2022 9.441185
ggplot(standard_dev_all_389_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 389 average temperatures for water years 1986-2021

MK & SS 389 (corrected)

sd_mk_389_ad <- mk.test(standard_dev_all_389_ad$sd_2)
print(sd_mk_389_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_389_ad$sd_2
## z = 1.0216, n = 32, p-value = 0.307
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   64.0000000 3802.6666667    0.1290323
sd_sens_389_ad <- sens.slope(standard_dev_all_389_ad$sd_2)
print(sd_sens_389_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_389_ad$sd_2
## z = 1.0216, n = 32, p-value = 0.307
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01140359  0.04896578
## sample estimates:
## Sen's slope 
##  0.01785583

Cumbres Trestle 431

Oyler -> Morrisey 6/8/2005

snotel_431 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "431")
#str(snotel_431) # check the date, usually a character.  

snotel_431$Date <- as.Date(snotel_431$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_431_clean <- snotel_431 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_431_clean <- snotel_431_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_431_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_431_clean <- snotel_431_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_min > -40) %>% 
  filter(temperature_min < 25)
ggplot(snotel_431_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
snotel_431_cull_count <- snotel_431_clean %>% 
  filter(temperature_min > -40) %>% 
  count(waterYear)

snotel_431_cull_count
## # A tibble: 36 x 2
## # Groups:   waterYear [36]
##    waterYear     n
##        <dbl> <int>
##  1      1987     1
##  2      1988   341
##  3      1989   364
##  4      1990   348
##  5      1991   341
##  6      1992   312
##  7      1993   356
##  8      1994   317
##  9      1995   364
## 10      1996   364
## # ... with 26 more rows
# filtering for too few observations in a year
snotel_431_cull_count_days <- snotel_431_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_431_cull_count_days
## # A tibble: 7 x 2
## # Groups:   waterYear [7]
##   waterYear     n
##       <dbl> <int>
## 1      1987     1
## 2      1988   341
## 3      1990   348
## 4      1991   341
## 5      1992   312
## 6      1994   317
## 7      2014   288

1987, 1988, 1990, 1991, 1992, 1994, 2014 need to be culled.

snotel_431_clean_culled <- snotel_431_clean %>% 
  filter(waterYear != "1987" & waterYear != "1988" & waterYear != "1990" & waterYear != "1991" & waterYear != "1992" & waterYear != "1994" & waterYear != "2014") #%>% 
  #filter(temperature_mean > -49)
ggplot(snotel_431_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_431_xts <- xts(snotel_431_clean_culled$temperature_mean, order.by = snotel_431_clean_culled$Date)

dygraph(temp_431_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
snotel_431_clean_culled <- snotel_431_clean_culled %>% 
  filter(temperature_mean < 20)

temp_431_xts <- xts(snotel_431_clean_culled$temperature_mean, order.by = snotel_431_clean_culled$Date)

dygraph(temp_431_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Cumbres Trestle 431

Oyler -> Morrisey 6/8/2005

snotel_431_adjusted <- snotel_431_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2005-06-08", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

431 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_431 <- snotel_431_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_431 <- yearly_wy_aver_431 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_431 <- daily_wy_aver_431 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_431$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_431 <-daily_wy_aver_431 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_431$date_temp <- signif(daily_wy_aver2_431$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_431, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

431 SD

standard_dev_431 <- daily_wy_aver_431 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_431 <- standard_dev_431 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_431 <- standard_dev_all_431 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_431 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1989 8.484011
1993 8.101376
1995 7.750446
1996 7.952021
1997 8.217323
1998 8.600343
1999 7.226028
2000 8.164156
2001 8.669911
2002 9.005161
2003 8.459072
2004 8.306872
2005 8.304679
2006 7.417238
2007 7.937180
2008 8.259566
2009 7.296965
2010 8.418225
2011 8.231369
2012 8.081428
2013 8.537331
2015 7.219164
2016 7.791453
2017 7.322790
2018 7.113402
2019 7.880447
2020 8.046531
2021 8.095523
2022 7.621146
ggplot(standard_dev_all_431, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 431 average temperatures for water years 2005-2021

MK & SS for 431 (non-corrected)

sd_mk_431 <- mk.test(standard_dev_all_431$sd_2)
print(sd_mk_431)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_431$sd_2
## z = -1.9321, n = 29, p-value = 0.05335
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -104.0000000 2842.0000000   -0.2561576
sd_sens_431 <- sens.slope(standard_dev_all_431$sd_2)
print(sd_sens_431)
## 
##  Sen's slope
## 
## data:  standard_dev_all_431$sd_2
## z = -1.9321, n = 29, p-value = 0.05335
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.041299369  0.001003281
## sample estimates:
## Sen's slope 
## -0.01890302

Corrected

431 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_431_ad <- snotel_431_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_431_ad <- yearly_wy_aver_431_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_431_ad <- daily_wy_aver_431_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_431_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_431_ad <-daily_wy_aver_431_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_431_ad$date_temp_ad <- signif(daily_wy_aver2_431_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_431_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

431 SS (corrected)

standard_dev_431_ad <- daily_wy_aver_431_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_431_ad <- standard_dev_431_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_431_ad <- standard_dev_all_431_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_431_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1989 7.914872
1993 7.518415
1995 7.173944
1996 7.365597
1997 7.634730
1998 7.970563
1999 6.682639
2000 7.537564
2001 8.039524
2002 8.344015
2003 7.813684
2004 7.737034
2005 7.599212
2006 7.418573
2007 7.938448
2008 8.260644
2009 7.298209
2010 8.419846
2011 8.232506
2012 8.082571
2013 8.538590
2015 7.220737
2016 7.793079
2017 7.324026
2018 7.114853
2019 7.881605
2020 8.047598
2021 8.096826
2022 7.622159
ggplot(standard_dev_all_431_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 431 average temperatures for water years 1986-2021

MK & SS 431 (corrected)

sd_mk_431_ad <- mk.test(standard_dev_all_431_ad$sd_2)
print(sd_mk_431_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_431_ad$sd_2
## z = 0.99418, n = 29, p-value = 0.3201
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   54.0000000 2842.0000000    0.1330049
sd_sens_431_ad <- sens.slope(standard_dev_all_431_ad$sd_2)
print(sd_sens_431_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_431_ad$sd_2
## z = 0.99418, n = 29, p-value = 0.3201
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01334615  0.03134200
## sample estimates:
## Sen's slope 
##  0.01029449

Idarado 538

Morrisey 7/12/2004

snotel_538 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "538")
#str(snotel_538) # check the date, usually a character.  

snotel_538$Date <- as.Date(snotel_538$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_538_clean <- snotel_538 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_538_clean <- snotel_538_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_538_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_538_clean <- snotel_538_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40)# %>% 
  #filter(temperature_min < 25)
ggplot(snotel_538_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
snotel_538_cull_count <- snotel_538_clean %>% 
  filter(temperature_min > -40) %>% 
  count(waterYear)

snotel_538_cull_count
## # A tibble: 36 x 2
## # Groups:   waterYear [36]
##    waterYear     n
##        <dbl> <int>
##  1      1987   355
##  2      1988   366
##  3      1989   364
##  4      1990   365
##  5      1991   337
##  6      1992   358
##  7      1993   346
##  8      1994   351
##  9      1995   364
## 10      1996   365
## # ... with 26 more rows
# filtering for too few observations in a year
snotel_538_cull_count_days <- snotel_538_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_538_cull_count_days
## # A tibble: 2 x 2
## # Groups:   waterYear [2]
##   waterYear     n
##       <dbl> <int>
## 1      1991   337
## 2      1993   346

1991 and 1993 need to be culled.

snotel_538_clean_culled <- snotel_538_clean %>% 
  filter(waterYear != "1991" & waterYear != "1993") # & waterYear != "1990" & waterYear != "1991" & waterYear != "1992" & waterYear != "1994" & waterYear != "2014") #%>% 
  #filter(temperature_mean > -49)
ggplot(snotel_538_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_538_xts <- xts(snotel_538_clean_culled$temperature_mean, order.by = snotel_538_clean_culled$Date)

dygraph(temp_538_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
snotel_538_clean_culled <- snotel_538_clean_culled %>% 
  filter(temperature_mean < 20)

temp_538_xts <- xts(snotel_538_clean_culled$temperature_mean, order.by = snotel_538_clean_culled$Date)

dygraph(temp_538_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Idarado 538

Morrisey 7/12/2004

snotel_538_adjusted <- snotel_538_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2004-07-12", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

538 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_538 <- snotel_538_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_538 <- yearly_wy_aver_538 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_538 <- daily_wy_aver_538 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_538$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_538 <-daily_wy_aver_538 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_538$date_temp <- signif(daily_wy_aver2_538$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_538, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

538 SD

standard_dev_538 <- daily_wy_aver_538 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_538 <- standard_dev_538 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_538 <- standard_dev_all_538 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_538 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 8.006364
1988 8.601864
1989 8.632529
1990 8.392986
1992 7.829013
1994 10.169815
1995 7.658097
1996 7.925011
1997 7.934668
1998 8.455869
1999 7.252626
2000 7.992580
2001 8.428201
2002 8.883209
2003 8.309823
2004 8.066471
2005 7.228205
2006 7.518896
2007 7.890984
2008 8.117130
2009 7.048527
2010 8.062733
2011 7.859864
2012 7.787906
2013 8.410412
2014 7.465408
2015 6.906228
2016 7.611126
2017 7.235848
2018 7.239780
2019 7.718304
2020 7.976082
2021 8.027478
2022 7.664387
ggplot(standard_dev_all_538, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 538 average temperatures for water years 2005-2021

MK & SS for 538 (non-corrected)

sd_mk_538 <- mk.test(standard_dev_all_538$sd_2)
print(sd_mk_538)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_538$sd_2
## z = -2.6684, n = 34, p-value = 0.007621
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -181.0000000 4550.3333333   -0.3226381
sd_sens_538 <- sens.slope(standard_dev_all_538$sd_2)
print(sd_sens_538)
## 
##  Sen's slope
## 
## data:  standard_dev_all_538$sd_2
## z = -2.6684, n = 34, p-value = 0.007621
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.042981556 -0.006896251
## sample estimates:
## Sen's slope 
## -0.02384088

Corrected

538 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_538_ad <- snotel_538_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_538_ad <- yearly_wy_aver_538_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_538_ad <- daily_wy_aver_538_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_538_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_538_ad <-daily_wy_aver_538_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_538_ad$date_temp_ad <- signif(daily_wy_aver2_538_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_538_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

538 SS (corrected)

standard_dev_538_ad <- daily_wy_aver_538_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_538_ad <- standard_dev_538_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_538_ad <- standard_dev_all_538_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_538_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 7.447584
1988 8.019241
1989 8.069873
1990 7.812557
1992 7.272544
1994 9.549437
1995 7.113641
1996 7.353811
1997 7.396458
1998 7.860736
1999 6.739299
2000 7.408528
2001 7.835100
2002 8.260161
2003 7.706375
2004 7.450714
2005 7.228187
2006 7.519259
2007 7.891254
2008 8.117349
2009 7.048722
2010 8.062844
2011 7.860215
2012 7.788573
2013 8.411054
2014 7.465688
2015 6.906947
2016 7.611742
2017 7.235817
2018 7.240024
2019 7.718445
2020 7.976612
2021 8.027693
2022 7.664991
ggplot(standard_dev_all_538_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 538 average temperatures for water years 1986-2021

MK & SS 538 (corrected)

sd_mk_538_ad <- mk.test(standard_dev_all_538_ad$sd_2)
print(sd_mk_538_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_538_ad$sd_2
## z = 0.14824, n = 34, p-value = 0.8821
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 1.100000e+01 4.550333e+03 1.960784e-02
sd_sens_538_ad <- sens.slope(standard_dev_all_538_ad$sd_2)
print(sd_sens_538_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_538_ad$sd_2
## z = 0.14824, n = 34, p-value = 0.8821
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01508416  0.01703536
## sample estimates:
##  Sen's slope 
## 0.0009105465

Lone Cone 589

Oyler -> Morrisey 6/22/2005

snotel_589 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "589")
#str(snotel_589) # check the date, usually a character.  

snotel_589$Date <- as.Date(snotel_589$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_589_clean <- snotel_589 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_589_clean <- snotel_589_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_589_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_589_clean <- snotel_589_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40)# %>% 
  #filter(temperature_min < 25)
ggplot(snotel_589_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
snotel_589_cull_count <- snotel_589_clean %>% 
  filter(temperature_min > -40) %>% 
  count(waterYear)

snotel_589_cull_count
## # A tibble: 37 x 2
## # Groups:   waterYear [37]
##    waterYear     n
##        <dbl> <int>
##  1      1986     1
##  2      1987   358
##  3      1988   366
##  4      1989   362
##  5      1990   365
##  6      1991   359
##  7      1992   366
##  8      1993   350
##  9      1994   337
## 10      1995   364
## # ... with 27 more rows
# filtering for too few observations in a year
snotel_589_cull_count_days <- snotel_589_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_589_cull_count_days
## # A tibble: 10 x 2
## # Groups:   waterYear [10]
##    waterYear     n
##        <dbl> <int>
##  1      1986     1
##  2      1994   337
##  3      1997   346
##  4      2006   317
##  5      2009   314
##  6      2010   332
##  7      2011   309
##  8      2012   322
##  9      2013   331
## 10      2014   335

1986, 1994, 1997, 2006, 2009, 2010, 2011, 2012, 2013, 2014 need to be culled.

snotel_589_clean_culled <- snotel_589_clean %>% 
  filter(waterYear != "1986" & waterYear != "1994" & waterYear != "1997" & waterYear != "2006" & waterYear != "2009" & waterYear != "2010" & waterYear != "2011" & waterYear != "2012" & waterYear != "2013" & waterYear != "2014") #%>% 
  #filter(temperature_mean > -49)
ggplot(snotel_589_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_589_xts <- xts(snotel_589_clean_culled$temperature_mean, order.by = snotel_589_clean_culled$Date)

dygraph(temp_589_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
snotel_589_clean_culled <- snotel_589_clean_culled %>% 
  filter(temperature_mean < 20)

temp_589_xts <- xts(snotel_589_clean_culled$temperature_mean, order.by = snotel_589_clean_culled$Date)

dygraph(temp_589_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Lone Cone 589

Oyler -> Morrisey 6/22/2005

snotel_589_adjusted <- snotel_589_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2005-06-22", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

589 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_589 <- snotel_589_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_589 <- yearly_wy_aver_589 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_589 <- daily_wy_aver_589 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_589$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_589 <-daily_wy_aver_589 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_589$date_temp <- signif(daily_wy_aver2_589$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_589, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

589 SD

standard_dev_589 <- daily_wy_aver_589 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_589 <- standard_dev_589 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_589 <- standard_dev_all_589 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_589 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 8.224944
1988 8.773281
1989 8.776693
1990 8.580289
1991 8.502502
1992 7.975611
1993 8.199081
1995 7.946318
1996 8.231469
1998 8.706109
1999 7.431979
2000 8.232988
2001 8.771085
2002 9.082003
2003 8.732769
2004 8.380419
2005 8.458900
2007 8.212611
2008 8.408220
2015 7.103666
2016 7.842413
2017 7.527960
2018 7.506884
2019 8.029199
2020 8.271099
2021 8.278625
2022 8.020852
ggplot(standard_dev_all_589, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 589 average temperatures for water years 2005-2021

MK & SS for 589 (non-corrected)

sd_mk_589 <- mk.test(standard_dev_all_589$sd_2)
print(sd_mk_589)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_589$sd_2
## z = -1.8345, n = 27, p-value = 0.06658
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  -89.0000000 2301.0000000   -0.2535613
sd_sens_589 <- sens.slope(standard_dev_all_589$sd_2)
print(sd_sens_589)
## 
##  Sen's slope
## 
## data:  standard_dev_all_589$sd_2
## z = -1.8345, n = 27, p-value = 0.06658
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.043902404  0.001230005
## sample estimates:
## Sen's slope 
## -0.02183399

Corrected

589 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_589_ad <- snotel_589_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_589_ad <- yearly_wy_aver_589_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_589_ad <- daily_wy_aver_589_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_589_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_589_ad <-daily_wy_aver_589_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_589_ad$date_temp_ad <- signif(daily_wy_aver2_589_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_589_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

589 SS (corrected)

standard_dev_589_ad <- daily_wy_aver_589_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_589_ad <- standard_dev_589_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_589_ad <- standard_dev_all_589_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_589_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 7.671688
1988 8.200568
1989 8.222864
1990 7.999469
1991 7.962227
1992 7.420415
1993 7.639311
1995 7.378994
1996 7.641722
1998 8.086353
1999 6.914950
2000 7.633161
2001 8.156564
2002 8.440080
2003 8.103839
2004 7.827141
2005 7.781158
2007 8.212720
2008 8.408841
2015 7.104293
2016 7.842763
2017 7.527954
2018 7.507329
2019 8.029667
2020 8.271310
2021 8.278726
2022 8.021229
ggplot(standard_dev_all_589_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 589 average temperatures for water years 1986-2021

MK & SS 589 (corrected)

sd_mk_589_ad <- mk.test(standard_dev_all_589_ad$sd_2)
print(sd_mk_589_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_589_ad$sd_2
## z = 0.6671, n = 27, p-value = 0.5047
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 3.300000e+01 2.301000e+03 9.401709e-02
sd_sens_589_ad <- sens.slope(standard_dev_all_589_ad$sd_2)
print(sd_sens_589_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_589_ad$sd_2
## z = 0.6671, n = 27, p-value = 0.5047
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01153608  0.02498422
## sample estimates:
## Sen's slope 
## 0.008250776

Middle Creek 624

Morrisey 6/26/2006

snotel_624 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "624")
#str(snotel_624) # check the date, usually a character.  

snotel_624$Date <- as.Date(snotel_624$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_624_clean <- snotel_624 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_624_clean <- snotel_624_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_624_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_624_clean <- snotel_624_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40) %>% 
  filter(temperature_mean < 25)
ggplot(snotel_624_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
snotel_624_cull_count <- snotel_624_clean %>% 
  filter(temperature_min > -40) %>% 
  count(waterYear)

snotel_624_cull_count
## # A tibble: 43 x 2
## # Groups:   waterYear [43]
##    waterYear     n
##        <dbl> <int>
##  1      1980    59
##  2      1981   354
##  3      1982     1
##  4      1983   312
##  5      1984   360
##  6      1985   362
##  7      1986   364
##  8      1987   360
##  9      1988   340
## 10      1989   363
## # ... with 33 more rows
# filtering for too few observations in a year
snotel_624_cull_count_days <- snotel_624_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_624_cull_count_days
## # A tibble: 9 x 2
## # Groups:   waterYear [9]
##   waterYear     n
##       <dbl> <int>
## 1      1980    59
## 2      1982     1
## 3      1983   312
## 4      1988   340
## 5      1991   338
## 6      1993   330
## 7      1994   330
## 8      2021   338
## 9      2022   348

1980, 1982, 1983, 1988, 1991, 1993, 1994, 2021, 2022 need to be culled.

snotel_624_clean_culled <- snotel_624_clean %>% 
  filter(waterYear != "1980" & waterYear != "1981" & waterYear != "1982" & waterYear != "1983" & waterYear != "1988" & waterYear != "1991" & waterYear != "1993" & waterYear != "1994" & waterYear != "2021" & waterYear != "2022")# & waterYear != "2014") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_624_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

ggplot(snotel_624_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_624_xts <- xts(snotel_624_clean_culled$temperature_mean, order.by = snotel_624_clean_culled$Date)

dygraph(temp_624_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
snotel_624_clean_culled <- snotel_624_clean_culled %>% 
  filter(temperature_mean < 20)

temp_624_xts <- xts(snotel_624_clean_culled$temperature_mean, order.by = snotel_624_clean_culled$Date)

dygraph(temp_624_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Middle Creek 624 Morrisey 6/26/2006

snotel_624_adjusted <- snotel_624_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2006-06-26", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

624 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_624 <- snotel_624_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_624 <- yearly_wy_aver_624 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_624 <- daily_wy_aver_624 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_624$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_624 <-daily_wy_aver_624 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_624$date_temp <- signif(daily_wy_aver2_624$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_624, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

624 SD

standard_dev_624 <- daily_wy_aver_624 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_624 <- standard_dev_624 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_624 <- standard_dev_all_624 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_624 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 8.239005
1985 8.207104
1986 7.191318
1987 7.871848
1989 8.454185
1990 8.200220
1992 6.438238
1995 7.565327
1996 7.835558
1997 7.996219
1998 8.415139
1999 7.174609
2000 7.839800
2001 8.355597
2002 8.492066
2003 8.150297
2004 7.912597
2005 7.915951
2006 7.936724
2007 7.640492
2008 8.021623
2009 7.171130
2010 8.259271
2011 7.932618
2012 7.791063
2013 8.259320
2014 7.372950
2015 6.806548
2016 7.548199
2017 7.108731
2018 7.036935
2019 7.675916
2020 7.741703
ggplot(standard_dev_all_624, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 624 average temperatures for water years 2005-2021

MK & SS for 624 (non-corrected)

sd_mk_624 <- mk.test(standard_dev_all_624$sd_2)
print(sd_mk_624)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_624$sd_2
## z = -1.9678, n = 33, p-value = 0.04909
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -128.0000000 4165.3333333   -0.2424242
sd_sens_624 <- sens.slope(standard_dev_all_624$sd_2)
print(sd_sens_624)
## 
##  Sen's slope
## 
## data:  standard_dev_all_624$sd_2
## z = -1.9678, n = 33, p-value = 0.04909
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.0333661824 -0.0003479011
## sample estimates:
## Sen's slope 
## -0.01662145

Corrected

624 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_624_ad <- snotel_624_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_624_ad <- yearly_wy_aver_624_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_624_ad <- daily_wy_aver_624_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_624_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_624_ad <-daily_wy_aver_624_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_624_ad$date_temp_ad <- signif(daily_wy_aver2_624_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_624_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

624 SS (corrected)

standard_dev_624_ad <- daily_wy_aver_624_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_624_ad <- standard_dev_624_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_624_ad <- standard_dev_all_624_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_624_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 7.738801
1985 7.713059
1986 6.725366
1987 7.374899
1989 7.942209
1990 7.687289
1992 6.080980
1995 7.080976
1996 7.336792
1997 7.519955
1998 7.885653
1999 6.718142
2000 7.329553
2001 7.833709
2002 7.944757
2003 7.617514
2004 7.442397
2005 7.390605
2006 7.304030
2007 7.639661
2008 8.021493
2009 7.170919
2010 8.259153
2011 7.932368
2012 7.790490
2013 8.259010
2014 7.372461
2015 6.806595
2016 7.547948
2017 7.107845
2018 7.036548
2019 7.675791
2020 7.741347
ggplot(standard_dev_all_624_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 624 average temperatures for water years 1986-2021

MK & SS 624 (corrected)

sd_mk_624_ad <- mk.test(standard_dev_all_624_ad$sd_2)
print(sd_mk_624_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_624_ad$sd_2
## z = 0.35637, n = 33, p-value = 0.7216
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 2.400000e+01 4.165333e+03 4.545455e-02
sd_sens_624_ad <- sens.slope(standard_dev_all_624_ad$sd_2)
print(sd_sens_624_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_624_ad$sd_2
## z = 0.35637, n = 33, p-value = 0.7216
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01304163  0.02223675
## sample estimates:
## Sen's slope 
## 0.003307866

Molas Lake 632

NSCE-Morrisey, Bias-Original 10/2/2003

snotel_632 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "632")
#str(snotel_632) # check the date, usually a character.  

snotel_632$Date <- as.Date(snotel_632$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_632_clean <- snotel_632 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_632_clean <- snotel_632_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_632_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_632_clean <- snotel_632_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  #filter(temperature_mean > -40) %>% 
  filter(temperature_mean < 25)
ggplot(snotel_632_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
snotel_632_cull_count <- snotel_632_clean %>% 
  filter(temperature_min > -40) %>% 
  count(waterYear)

snotel_632_cull_count
## # A tibble: 37 x 2
## # Groups:   waterYear [37]
##    waterYear     n
##        <dbl> <int>
##  1      1986    54
##  2      1987   355
##  3      1988   362
##  4      1989   364
##  5      1990   365
##  6      1991   365
##  7      1992   366
##  8      1993   303
##  9      1994   294
## 10      1995   361
## # ... with 27 more rows
# filtering for too few observations in a year
snotel_632_cull_count_days <- snotel_632_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_632_cull_count_days
## # A tibble: 6 x 2
## # Groups:   waterYear [6]
##   waterYear     n
##       <dbl> <int>
## 1      1986    54
## 2      1993   303
## 3      1994   295
## 4      2005   344
## 5      2006   304
## 6      2013   248

1986, 1993, 1994, 2005, 2006, 2013 need to be culled.

snotel_632_clean_culled <- snotel_632_clean %>% 
  filter(waterYear != "1986" & waterYear != "1993" & waterYear != "1994" & waterYear != "2005" & waterYear != "2006" & waterYear != "2013")# & waterYear != "1993" & waterYear != "1994" & waterYear != "2021" & waterYear != "2022") & waterYear != "2014") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_632_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

ggplot(snotel_632_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_632_xts <- xts(snotel_632_clean_culled$temperature_mean, order.by = snotel_632_clean_culled$Date)

dygraph(temp_632_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
snotel_632_clean_culled <- snotel_632_clean_culled %>% 
  filter(temperature_mean < 20)

temp_632_xts <- xts(snotel_632_clean_culled$temperature_mean, order.by = snotel_632_clean_culled$Date)

dygraph(temp_632_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Molas Lake 632

NSCE-Morrisey, Bias-Original 10/2/2003

snotel_632_adjusted <- snotel_632_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2003-10-02", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

632 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_632 <- snotel_632_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_632 <- yearly_wy_aver_632 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_632 <- daily_wy_aver_632 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_632$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_632 <-daily_wy_aver_632 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_632$date_temp <- signif(daily_wy_aver2_632$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_632, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

632 SD

standard_dev_632 <- daily_wy_aver_632 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_632 <- standard_dev_632 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_632 <- standard_dev_all_632 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_632 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 7.613234
1988 8.191410
1989 8.166070
1990 7.646926
1991 7.969691
1992 7.156808
1995 7.489888
1996 7.655841
1997 7.698012
1998 8.157837
1999 6.885991
2000 7.657554
2001 6.788385
2002 8.254532
2003 7.907020
2004 7.308264
2007 7.607147
2008 7.964946
2009 6.957483
2010 7.983696
2011 7.701623
2012 7.490054
2014 7.236529
2015 6.581104
2016 7.409875
2017 7.001441
2018 6.814621
2019 7.582515
2020 7.766470
2021 7.792571
2022 7.357711
ggplot(standard_dev_all_632, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 632 average temperatures for water years 2005-2021

MK & SS for 632 (non-corrected)

sd_mk_632 <- mk.test(standard_dev_all_632$sd_2)
print(sd_mk_632)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_632$sd_2
## z = -1.6317, n = 31, p-value = 0.1028
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  -97.0000000 3461.6666667   -0.2086022
sd_sens_632 <- sens.slope(standard_dev_all_632$sd_2)
print(sd_sens_632)
## 
##  Sen's slope
## 
## data:  standard_dev_all_632$sd_2
## z = -1.6317, n = 31, p-value = 0.1028
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.034372524  0.003422909
## sample estimates:
## Sen's slope 
## -0.01583493

Corrected

632 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_632_ad <- snotel_632_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_632_ad <- yearly_wy_aver_632_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_632_ad <- daily_wy_aver_632_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_632_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_632_ad <-daily_wy_aver_632_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_632_ad$date_temp_ad <- signif(daily_wy_aver2_632_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_632_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

632 SS (corrected)

standard_dev_632_ad <- daily_wy_aver_632_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_632_ad <- standard_dev_632_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_632_ad <- standard_dev_all_632_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_632_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 7.106416
1988 7.678134
1989 7.662133
1990 7.141729
1991 7.492618
1992 6.671512
1995 6.957618
1996 7.110780
1997 7.178235
1998 7.586041
1999 6.397699
2000 7.103150
2001 6.358660
2002 7.672831
2003 7.338249
2004 7.309361
2007 7.606845
2008 7.965266
2009 6.957636
2010 7.984148
2011 7.701741
2012 7.490245
2014 7.236308
2015 6.581364
2016 7.409996
2017 7.000722
2018 6.815246
2019 7.583475
2020 7.766635
2021 7.793272
2022 7.358050
ggplot(standard_dev_all_632_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 632 average temperatures for water years 1986-2021

MK & SS 632 (corrected)

sd_mk_632_ad <- mk.test(standard_dev_all_632_ad$sd_2)
print(sd_mk_632_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_632_ad$sd_2
## z = 0.9518, n = 31, p-value = 0.3412
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   57.0000000 3461.6666667    0.1225806
sd_sens_632_ad <- sens.slope(standard_dev_all_632_ad$sd_2)
print(sd_sens_632_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_632_ad$sd_2
## z = 0.9518, n = 31, p-value = 0.3412
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01086011  0.02505935
## sample estimates:
## Sen's slope 
## 0.008387789

Red Mountain Pass 713

Morrisey 8/18/2004

snotel_713 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "713")
#str(snotel_713) # check the date, usually a character.  

snotel_713$Date <- as.Date(snotel_713$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_713_clean <- snotel_713 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_713_clean <- snotel_713_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_713_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_713_clean <- snotel_713_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40) %>% 
  filter(temperature_mean < 25)
ggplot(snotel_713_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
snotel_713_cull_count <- snotel_713_clean %>% 
  filter(temperature_min > -40) %>% 
  count(waterYear)

snotel_713_cull_count
## # A tibble: 42 x 2
## # Groups:   waterYear [42]
##    waterYear     n
##        <dbl> <int>
##  1      1981   189
##  2      1982   287
##  3      1983   319
##  4      1984   360
##  5      1985   362
##  6      1986   364
##  7      1987   263
##  8      1988   282
##  9      1989   228
## 10      1990   365
## # ... with 32 more rows
# filtering for too few observations in a year
snotel_713_cull_count_days <- snotel_713_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_713_cull_count_days
## # A tibble: 6 x 2
## # Groups:   waterYear [6]
##   waterYear     n
##       <dbl> <int>
## 1      1981   189
## 2      1982   287
## 3      1983   319
## 4      1987   315
## 5      1989   263
## 6      1994   308

1981, 1982, 1983, 1987, 1989, 1994 need to be culled. 1984 & 1985 are too low.

snotel_713_clean_culled <- snotel_713_clean %>% 
  filter(waterYear > "1985" & waterYear != "1987" & waterYear != "1989" & waterYear != "1994")# & waterYear != "1993" & waterYear != "1994" & waterYear != "2021" & waterYear != "2022") & waterYear != "2014") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_713_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

ggplot(snotel_713_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_713_xts <- xts(snotel_713_clean_culled$temperature_mean, order.by = snotel_713_clean_culled$Date)

dygraph(temp_713_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
snotel_713_clean_culled <- snotel_713_clean_culled %>% 
  filter(temperature_mean < 19.3)

temp_713_xts <- xts(snotel_713_clean_culled$temperature_mean, order.by = snotel_713_clean_culled$Date)

dygraph(temp_713_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Red Mountain Pass 713

Morrisey 8/18/2004

snotel_713_adjusted <- snotel_713_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2004-08-18", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

713 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_713 <- snotel_713_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_713 <- yearly_wy_aver_713 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_713 <- daily_wy_aver_713 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_713$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_713 <-daily_wy_aver_713 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_713$date_temp <- signif(daily_wy_aver2_713$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_713, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

713 SD

standard_dev_713 <- daily_wy_aver_713 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_713 <- standard_dev_713 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_713 <- standard_dev_all_713 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_713 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1986 7.022683
1988 9.344095
1990 7.761743
1991 7.768555
1992 7.212083
1993 7.716280
1995 7.336781
1996 7.731429
1997 7.947757
1998 8.198019
1999 6.921528
2000 7.680921
2001 8.089853
2002 8.395619
2003 7.952564
2004 7.777488
2005 7.151965
2006 7.312223
2007 7.636775
2008 7.908986
2009 6.940732
2010 7.877650
2011 7.730288
2012 7.531887
2013 8.195738
2014 7.284701
2015 6.638343
2016 7.380700
2017 7.097935
2018 7.002862
2019 7.626598
2020 7.810018
2021 7.817233
2022 7.339213
ggplot(standard_dev_all_713, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 713 average temperatures for water years 2005-2021

MK & SS for 713 (non-corrected)

sd_mk_713 <- mk.test(standard_dev_all_713$sd_2)
print(sd_mk_713)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_713$sd_2
## z = -1.1563, n = 34, p-value = 0.2476
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##          S       varS        tau 
##  -79.00000 4550.33333   -0.14082
sd_sens_713 <- sens.slope(standard_dev_all_713$sd_2)
print(sd_sens_713)
## 
##  Sen's slope
## 
## data:  standard_dev_all_713$sd_2
## z = -1.1563, n = 34, p-value = 0.2476
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.028289222  0.006454843
## sample estimates:
##  Sen's slope 
## -0.008980292

Corrected

713 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_713_ad <- snotel_713_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_713_ad <- yearly_wy_aver_713_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_713_ad <- daily_wy_aver_713_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_713_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_713_ad <-daily_wy_aver_713_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_713_ad$date_temp_ad <- signif(daily_wy_aver2_713_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_713_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

713 SS (corrected)

standard_dev_713_ad <- daily_wy_aver_713_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_713_ad <- standard_dev_713_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_713_ad <- standard_dev_all_713_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_713_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1986 6.564707
1988 8.872885
1990 7.267850
1991 7.311859
1992 6.739882
1993 7.221869
1995 6.846024
1996 7.223659
1997 7.456352
1998 7.661532
1999 6.466454
2000 7.159552
2001 7.557097
2002 7.840372
2003 7.411227
2004 7.238021
2005 7.151835
2006 7.312105
2007 7.636625
2008 7.909548
2009 6.940894
2010 7.878182
2011 7.730889
2012 7.532385
2013 8.196125
2014 7.284711
2015 6.639012
2016 7.380681
2017 7.097919
2018 7.003100
2019 7.627224
2020 7.810395
2021 7.817455
2022 7.339782
ggplot(standard_dev_all_713_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 713 average temperatures for water years 1986-2021

MK & SS 713 (corrected)

sd_mk_713_ad <- mk.test(standard_dev_all_713_ad$sd_2)
print(sd_mk_713_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_713_ad$sd_2
## z = 1.4231, n = 34, p-value = 0.1547
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   97.0000000 4550.3333333    0.1729055
sd_sens_713_ad <- sens.slope(standard_dev_all_713_ad$sd_2)
print(sd_sens_713_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_713_ad$sd_2
## z = 1.4231, n = 34, p-value = 0.1547
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.005520634  0.030221234
## sample estimates:
## Sen's slope 
##  0.01194811

Scotch Creek 739

Morrisey 6/15/2004

snotel_739 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "739")
#str(snotel_739) # check the date, usually a character.  

snotel_739$Date <- as.Date(snotel_739$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_739_clean <- snotel_739 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_739_clean <- snotel_739_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_739_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_739_clean <- snotel_739_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40) %>% 
  filter(temperature_mean < 25)
ggplot(snotel_739_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
snotel_739_cull_count <- snotel_739_clean %>% 
  filter(temperature_min > -40) %>% 
  count(waterYear)

snotel_739_cull_count
## # A tibble: 36 x 2
## # Groups:   waterYear [36]
##    waterYear     n
##        <dbl> <int>
##  1      1987   345
##  2      1988   273
##  3      1989   363
##  4      1990   365
##  5      1991   365
##  6      1992   365
##  7      1993   342
##  8      1994   345
##  9      1995   362
## 10      1996   365
## # ... with 26 more rows
# filtering for too few observations in a year
snotel_739_cull_count_days <- snotel_739_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_739_cull_count_days
## # A tibble: 6 x 2
## # Groups:   waterYear [6]
##   waterYear     n
##       <dbl> <int>
## 1      1987   345
## 2      1988   273
## 3      1993   343
## 4      1994   345
## 5      2016   320
## 6      2017   249

1987, 1988, 1993, 1994, 2016, 2017 need to be culled.

snotel_739_clean_culled <- snotel_739_clean %>% 
  filter(waterYear != "1987" & waterYear != "1988" & waterYear != "1993" & waterYear != "1994" & waterYear != "2016" & waterYear != "2017")# & waterYear != "1993" & waterYear != "1994" & waterYear != "2021" & waterYear != "2022") & waterYear != "2014") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_739_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

ggplot(snotel_739_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_739_xts <- xts(snotel_739_clean_culled$temperature_mean, order.by = snotel_739_clean_culled$Date)

dygraph(temp_739_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
snotel_739_clean_culled <- snotel_739_clean_culled %>% 
  filter(temperature_mean < 19.3)

temp_739_xts <- xts(snotel_739_clean_culled$temperature_mean, order.by = snotel_739_clean_culled$Date)

dygraph(temp_739_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Scotch Creek 739

Morrisey 6/15/2004

snotel_739_adjusted <- snotel_739_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2004-06-15", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

739 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_739 <- snotel_739_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_739 <- yearly_wy_aver_739 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_739 <- daily_wy_aver_739 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_739$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_739 <-daily_wy_aver_739 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_739$date_temp <- signif(daily_wy_aver2_739$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_739, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

739 SD

standard_dev_739 <- daily_wy_aver_739 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_739 <- standard_dev_739 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_739 <- standard_dev_all_739 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_739 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1989 7.327097
1990 6.021389
1991 5.989983
1992 4.863013
1995 7.600315
1996 7.840086
1997 7.952598
1998 8.393336
1999 7.124567
2000 7.697346
2001 8.108824
2002 8.727256
2003 7.977395
2004 7.894210
2005 7.144352
2006 7.250093
2007 7.615681
2008 7.916198
2009 6.788880
2010 7.880288
2011 7.703849
2012 7.458561
2013 8.171355
2014 7.223604
2015 6.545241
2018 6.908918
2019 7.558347
2020 7.673725
2021 7.658883
2022 7.427996
ggplot(standard_dev_all_739, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 739 average temperatures for water years 2005-2021

MK & SS for 739 (non-corrected)

sd_mk_739 <- mk.test(standard_dev_all_739$sd_2)
print(sd_mk_739)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_739$sd_2
## z = 0, n = 30, p-value = 1
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 1.000000e+00 3.141667e+03 2.298851e-03
sd_sens_739 <- sens.slope(standard_dev_all_739$sd_2)
print(sd_sens_739)
## 
##  Sen's slope
## 
## data:  standard_dev_all_739$sd_2
## z = 0, n = 30, p-value = 1
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.02559406  0.03474477
## sample estimates:
##  Sen's slope 
## 0.0005911373

Corrected

739 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_739_ad <- snotel_739_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_739_ad <- yearly_wy_aver_739_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_739_ad <- daily_wy_aver_739_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_739_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_739_ad <-daily_wy_aver_739_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_739_ad$date_temp_ad <- signif(daily_wy_aver2_739_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_739_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

739 SS (corrected)

standard_dev_739_ad <- daily_wy_aver_739_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_739_ad <- standard_dev_739_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_739_ad <- standard_dev_all_739_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_739_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1989 6.882254
1990 5.641380
1991 5.682802
1992 4.573743
1995 7.029026
1996 7.246834
1997 7.374205
1998 7.773434
1999 6.583515
2000 7.110733
2001 7.493196
2002 8.076937
2003 7.362576
2004 7.247292
2005 7.144310
2006 7.250207
2007 7.615515
2008 7.916383
2009 6.789140
2010 7.880400
2011 7.703781
2012 7.458914
2013 8.171585
2014 7.223615
2015 6.545250
2018 6.908860
2019 7.558520
2020 7.674026
2021 7.658936
2022 7.428362
ggplot(standard_dev_all_739_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 739 average temperatures for water years 1986-2021

MK & SS 739 (corrected)

sd_mk_739_ad <- mk.test(standard_dev_all_739_ad$sd_2)
print(sd_mk_739_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_739_ad$sd_2
## z = 2.3907, n = 30, p-value = 0.01682
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  135.0000000 3141.6666667    0.3103448
sd_sens_739_ad <- sens.slope(standard_dev_all_739_ad$sd_2)
print(sd_sens_739_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_739_ad$sd_2
## z = 2.3907, n = 30, p-value = 0.01682
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  0.00806300 0.05860595
## sample estimates:
## Sen's slope 
##  0.02802846

Slumgullion 762

Morrisey 6/26/2006

snotel_762 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "762")
#str(snotel_762) # check the date, usually a character.  

snotel_762$Date <- as.Date(snotel_762$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_762_clean <- snotel_762 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_762_clean <- snotel_762_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_762_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_762_clean <- snotel_762_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40) %>% 
  filter(temperature_mean < 25)
ggplot(snotel_762_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
#snotel_762_cull_count <- snotel_762_clean %>% 
#  filter(temperature_min > -40) %>% 
#  count(waterYear)

#snotel_762_cull_count

# filtering for too few observations in a year
snotel_762_cull_count_days <- snotel_762_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_762_cull_count_days
## # A tibble: 5 x 2
## # Groups:   waterYear [5]
##   waterYear     n
##       <dbl> <int>
## 1      1980    43
## 2      1982     1
## 3      1983   319
## 4      1994   331
## 5      2003   337

1980, 1982, 1994, 2003 need to be culled. 1995 maximums are less than the minimums.

snotel_762_clean_culled <- snotel_762_clean %>% 
  filter(waterYear != "1980" & waterYear != "1981" & waterYear != "1982" & waterYear != "1983" & waterYear != "1994" & waterYear != "1995" & waterYear != "2003")# & waterYear != "2017" & waterYear != "1993" & waterYear != "1994" & waterYear != "2021" & waterYear != "2022") & waterYear != "2014") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_762_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

ggplot(snotel_762_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_762_xts <- xts(snotel_762_clean_culled$temperature_mean, order.by = snotel_762_clean_culled$Date)

dygraph(temp_762_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
#snotel_762_clean_culled <- snotel_762_clean_culled %>% 
#  filter(temperature_mean < 19.3)

temp_762_xts <- xts(snotel_762_clean_culled$temperature_mean, order.by = snotel_762_clean_culled$Date)

dygraph(temp_762_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Slumgullion 762

Morrisey 6/26/2006

snotel_762_adjusted <- snotel_762_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2006-06-26", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

762 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_762 <- snotel_762_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_762 <- yearly_wy_aver_762 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_762 <- daily_wy_aver_762 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_762$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_762 <-daily_wy_aver_762 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_762$date_temp <- signif(daily_wy_aver2_762$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_762, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

762 SD

standard_dev_762 <- daily_wy_aver_762 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_762 <- standard_dev_762 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_762 <- standard_dev_all_762 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_762 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 8.474614
1985 8.395416
1986 7.224771
1987 8.036169
1988 8.470217
1989 8.556810
1990 8.104917
1991 8.129177
1992 7.378858
1993 7.881362
1996 8.151725
1997 8.002402
1998 8.492764
1999 7.140044
2000 7.859908
2001 8.456394
2002 8.777834
2004 8.172096
2005 8.157803
2006 8.479381
2007 7.832136
2008 8.378519
2009 7.328609
2010 8.218219
2011 8.186856
2012 7.846460
2013 8.493891
2014 7.772861
2015 7.080609
2016 7.803357
2017 7.532567
2018 7.423620
2019 8.010211
2020 8.142119
2021 8.245368
2022 7.614681
ggplot(standard_dev_all_762, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 762 average temperatures for water years 2005-2021

MK & SS for 762 (non-corrected)

sd_mk_762 <- mk.test(standard_dev_all_762$sd_2)
print(sd_mk_762)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_762$sd_2
## z = -1.4574, n = 36, p-value = 0.145
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -108.0000000 5390.0000000   -0.1714286
sd_sens_762 <- sens.slope(standard_dev_all_762$sd_2)
print(sd_sens_762)
## 
##  Sen's slope
## 
## data:  standard_dev_all_762$sd_2
## z = -1.4574, n = 36, p-value = 0.145
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.024958596  0.002910121
## sample estimates:
## Sen's slope 
## -0.01112266

Corrected

762 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_762_ad <- snotel_762_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_762_ad <- yearly_wy_aver_762_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_762_ad <- daily_wy_aver_762_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_762_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_762_ad <-daily_wy_aver_762_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_762_ad$date_temp_ad <- signif(daily_wy_aver2_762_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_762_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

762 SS (corrected)

standard_dev_762_ad <- daily_wy_aver_762_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_762_ad <- standard_dev_762_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_762_ad <- standard_dev_all_762_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_762_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 7.970791
1985 7.888311
1986 6.764665
1987 7.528970
1988 7.965877
1989 8.053090
1990 7.604786
1991 7.665906
1992 6.904408
1993 7.400370
1996 7.629529
1997 7.511155
1998 7.956453
1999 6.688385
2000 7.353371
2001 7.919104
2002 8.206929
2004 7.685091
2005 7.619040
2006 7.823927
2007 7.831876
2008 8.379138
2009 7.329078
2010 8.218734
2011 8.187217
2012 7.846411
2013 8.493922
2014 7.773398
2015 7.081459
2016 7.804132
2017 7.532344
2018 7.424275
2019 8.010639
2020 8.142599
2021 8.245690
2022 7.615524
ggplot(standard_dev_all_762_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 762 average temperatures for water years 1986-2021

MK & SS 762 (corrected)

sd_mk_762_ad <- mk.test(standard_dev_all_762_ad$sd_2)
print(sd_mk_762_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_762_ad$sd_2
## z = 1.2123, n = 36, p-value = 0.2254
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   90.0000000 5390.0000000    0.1428571
sd_sens_762_ad <- sens.slope(standard_dev_all_762_ad$sd_2)
print(sd_sens_762_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_762_ad$sd_2
## z = 1.2123, n = 36, p-value = 0.2254
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.00495755  0.02329061
## sample estimates:
## Sen's slope 
## 0.009103679

Spud Mountain 780

Morrisey 6/28/2004

snotel_780 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "780")
#str(snotel_780) # check the date, usually a character.  

snotel_780$Date <- as.Date(snotel_780$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_780_clean <- snotel_780 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_780_clean <- snotel_780_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_780_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_780_clean <- snotel_780_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40) %>% 
  filter(temperature_mean < 25)
ggplot(snotel_780_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
#snotel_780_cull_count <- snotel_780_clean %>% 
#  filter(temperature_min > -40) %>% 
#  count(waterYear)

#snotel_780_cull_count

# filtering for too few observations in a year
snotel_780_cull_count_days <- snotel_780_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_780_cull_count_days
## # A tibble: 8 x 2
## # Groups:   waterYear [8]
##   waterYear     n
##       <dbl> <int>
## 1      1987   315
## 2      1994   336
## 3      1998   309
## 4      2009   342
## 5      2011   337
## 6      2013   336
## 7      2016   331
## 8      2017   324

1987, 1994, 1998, 2009, 2011, 2013, 2016, 2017 need to be culled.

snotel_780_clean_culled <- snotel_780_clean %>% 
  filter(waterYear != "1987" & waterYear != "1994" & waterYear != "1998" & waterYear != "2009" & waterYear != "2011" & waterYear != "2013" & waterYear != "2016" & waterYear != "2017")# & waterYear != "1993" & waterYear != "1994" & waterYear != "2021" & waterYear != "2022") & waterYear != "2014") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_780_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

ggplot(snotel_780_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_780_xts <- xts(snotel_780_clean_culled$temperature_mean, order.by = snotel_780_clean_culled$Date)

dygraph(temp_780_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
#snotel_780_clean_culled <- snotel_780_clean_culled %>% 
#  filter(temperature_mean < 19.3)

temp_780_xts <- xts(snotel_780_clean_culled$temperature_mean, order.by = snotel_780_clean_culled$Date)

dygraph(temp_780_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Spud Mountain 780

Morrisey 6/28/2004

snotel_780_adjusted <- snotel_780_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2004-06-28", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

780 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_780 <- snotel_780_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_780 <- yearly_wy_aver_780 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_780 <- daily_wy_aver_780 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_780$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_780 <-daily_wy_aver_780 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_780$date_temp <- signif(daily_wy_aver2_780$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_780, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

780 SD

standard_dev_780 <- daily_wy_aver_780 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_780 <- standard_dev_780 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_780 <- standard_dev_all_780 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_780 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1988 8.083874
1989 8.173292
1990 7.631378
1991 7.796215
1992 7.011676
1993 7.893616
1995 7.478626
1996 7.773665
1997 7.856295
1999 6.828980
2000 7.722333
2001 8.147987
2002 8.545717
2003 7.955225
2004 7.944675
2005 7.146399
2006 7.088438
2007 7.424508
2008 7.859459
2010 7.733705
2012 7.461301
2014 7.143839
2015 6.533563
2018 6.864656
2019 7.580081
2020 7.727076
2021 7.724978
2022 7.146411
ggplot(standard_dev_all_780, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 780 average temperatures for water years 2005-2021

MK & SS for 780 (non-corrected)

sd_mk_780 <- mk.test(standard_dev_all_780$sd_2)
print(sd_mk_780)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_780$sd_2
## z = -1.9954, n = 28, p-value = 0.046
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -102.0000000 2562.0000000   -0.2698413
sd_sens_780 <- sens.slope(standard_dev_all_780$sd_2)
print(sd_sens_780)
## 
##  Sen's slope
## 
## data:  standard_dev_all_780$sd_2
## z = -1.9954, n = 28, p-value = 0.046
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.0394954014 -0.0004267439
## sample estimates:
## Sen's slope 
## -0.01908281

Corrected

780 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_780_ad <- snotel_780_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_780_ad <- yearly_wy_aver_780_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_780_ad <- daily_wy_aver_780_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_780_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_780_ad <-daily_wy_aver_780_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_780_ad$date_temp_ad <- signif(daily_wy_aver2_780_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_780_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

780 SS (corrected)

standard_dev_780_ad <- daily_wy_aver_780_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_780_ad <- standard_dev_780_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_780_ad <- standard_dev_all_780_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_780_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1988 7.530892
1989 7.629513
1990 7.085085
1991 7.275264
1992 6.494536
1993 7.336765
1995 6.928211
1996 7.209777
1997 7.307706
1999 6.331029
2000 7.141268
2001 7.555861
2002 7.925084
2003 7.356618
2004 7.299749
2005 7.146735
2006 7.089360
2007 7.425654
2008 7.860535
2010 7.734756
2012 7.462985
2014 7.144698
2015 6.534943
2018 6.866003
2019 7.581136
2020 7.728271
2021 7.725902
2022 7.147773
ggplot(standard_dev_all_780_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 780 average temperatures for water years 1986-2021

MK & SS 780 (corrected)

sd_mk_780_ad <- mk.test(standard_dev_all_780_ad$sd_2)
print(sd_mk_780_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_780_ad$sd_2
## z = 0.84953, n = 28, p-value = 0.3956
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   44.0000000 2562.0000000    0.1164021
sd_sens_780_ad <- sens.slope(standard_dev_all_780_ad$sd_2)
print(sd_sens_780_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_780_ad$sd_2
## z = 0.84953, n = 28, p-value = 0.3956
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01274116  0.02720599
## sample estimates:
## Sen's slope 
## 0.007453894

Stump Lakes 797

Morrisey 7/22/2005

snotel_797 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "797")
#str(snotel_797) # check the date, usually a character.  

snotel_797$Date <- as.Date(snotel_797$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_797_clean <- snotel_797 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_797_clean <- snotel_797_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_797_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_797_clean <- snotel_797_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40) %>% 
  filter(temperature_mean < 40)
ggplot(snotel_797_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
#snotel_797_cull_count <- snotel_797_clean %>% 
#  filter(temperature_min > -40) %>% 
#  count(waterYear)

#snotel_797_cull_count

# filtering for too few observations in a year
snotel_797_cull_count_days <- snotel_797_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_797_cull_count_days
## # A tibble: 5 x 2
## # Groups:   waterYear [5]
##   waterYear     n
##       <dbl> <int>
## 1      1987   267
## 2      1993   320
## 3      1994   274
## 4      2005   329
## 5      2006   327

1987, 1993, 1994, 2005, 2006 need to be culled.

snotel_797_clean_culled <- snotel_797_clean %>% 
  filter(waterYear != "1994" & waterYear != "2005" & waterYear != "2006" & waterYear != "1993" & waterYear != "1987")# & waterYear != "2021" & waterYear != "2022") & waterYear != "2014") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_797_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

ggplot(snotel_797_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_797_xts <- xts(snotel_797_clean_culled$temperature_mean, order.by = snotel_797_clean_culled$Date)

dygraph(temp_797_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
#snotel_797_clean_culled <- snotel_797_clean_culled %>% 
#  filter(temperature_mean < 19.3)

temp_797_xts <- xts(snotel_797_clean_culled$temperature_mean, order.by = snotel_797_clean_culled$Date)

dygraph(temp_797_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Stump Lakes 797

Morrisey 7/22/2005

snotel_797_adjusted <- snotel_797_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2005-07-22", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

797 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_797 <- snotel_797_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_797 <- yearly_wy_aver_797 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_797 <- daily_wy_aver_797 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_797$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_797 <-daily_wy_aver_797 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_797$date_temp <- signif(daily_wy_aver2_797$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_797, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

797 SD

standard_dev_797 <- daily_wy_aver_797 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_797 <- standard_dev_797 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_797 <- standard_dev_all_797 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_797 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1988 8.133643
1989 10.300128
1990 7.808751
1991 9.471092
1992 7.312520
1995 7.694636
1996 7.953496
1997 8.069085
1998 8.367447
1999 7.015052
2000 7.935172
2001 8.545423
2002 9.270169
2003 8.084230
2004 9.099232
2007 7.603452
2008 8.008843
2009 7.184718
2010 7.913867
2011 7.886632
2012 7.696409
2013 8.096842
2014 7.327713
2015 6.737528
2016 7.565177
2017 7.316786
2018 7.076498
2019 7.800277
2020 7.872608
2021 7.863889
2022 7.341735
ggplot(standard_dev_all_797, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 797 average temperatures for water years 2005-2021

MK & SS for 797 (non-corrected)

sd_mk_797 <- mk.test(standard_dev_all_797$sd_2)
print(sd_mk_797)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_797$sd_2
## z = -2.5495, n = 31, p-value = 0.01079
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -151.0000000 3461.6666667   -0.3247312
sd_sens_797 <- sens.slope(standard_dev_all_797$sd_2)
print(sd_sens_797)
## 
##  Sen's slope
## 
## data:  standard_dev_all_797$sd_2
## z = -2.5495, n = 31, p-value = 0.01079
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.055919889 -0.005393295
## sample estimates:
## Sen's slope 
## -0.02911914

Corrected

797 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_797_ad <- snotel_797_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_797_ad <- yearly_wy_aver_797_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_797_ad <- daily_wy_aver_797_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_797_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_797_ad <-daily_wy_aver_797_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_797_ad$date_temp_ad <- signif(daily_wy_aver2_797_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_797_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

797 SS (corrected)

standard_dev_797_ad <- daily_wy_aver_797_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_797_ad <- standard_dev_797_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_797_ad <- standard_dev_all_797_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_797_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1988 7.629765
1989 9.567031
1990 7.287673
1991 8.816683
1992 6.729102
1995 7.165683
1996 7.422437
1997 7.545973
1998 7.813938
1999 6.542177
2000 7.377392
2001 7.978790
2002 8.614161
2003 7.487087
2004 8.402698
2007 7.603448
2008 8.008551
2009 7.184680
2010 7.914300
2011 7.886697
2012 7.696830
2013 8.096444
2014 7.327715
2015 6.737458
2016 7.565155
2017 7.316299
2018 7.076457
2019 7.800532
2020 7.873011
2021 7.863666
2022 7.341402
ggplot(standard_dev_all_797_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 797 average temperatures for water years 1986-2021

MK & SS 797 (corrected)

sd_mk_797_ad <- mk.test(standard_dev_all_797_ad$sd_2)
print(sd_mk_797_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_797_ad$sd_2
## z = -0.30594, n = 31, p-value = 0.7597
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##  -19.00000000 3461.66666667   -0.04086022
sd_sens_797_ad <- sens.slope(standard_dev_all_797_ad$sd_2)
print(sd_sens_797_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_797_ad$sd_2
## z = -0.30594, n = 31, p-value = 0.7597
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.03291859  0.02011237
## sample estimates:
## Sen's slope 
## -0.00412891

Upper Rio Grande 839

Oyler -> Morrisey 5/26/2004

snotel_839 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "839")
#str(snotel_839) # check the date, usually a character.  

snotel_839$Date <- as.Date(snotel_839$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_839_clean <- snotel_839 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_839_clean <- snotel_839_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_839_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_839_clean <- snotel_839_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40) %>% 
  filter(temperature_mean < 25)
ggplot(snotel_839_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
#snotel_839_cull_count <- snotel_839_clean %>% 
#  filter(temperature_min > -40) %>% 
#  count(waterYear)

#snotel_839_cull_count

# filtering for too few observations in a year
snotel_839_cull_count_days <- snotel_839_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_839_cull_count_days
## # A tibble: 5 x 2
## # Groups:   waterYear [5]
##   waterYear     n
##       <dbl> <int>
## 1      1994   338
## 2      2004   330
## 3      2012   339
## 4      2014   313
## 5      2016   347

1994, 2004, 2012, 2014, 2016 need to be culled.

snotel_839_clean_culled <- snotel_839_clean %>% 
  filter(waterYear != "1994" & waterYear != "2004" & waterYear != "2012" & waterYear != "2014" & waterYear != "2016")# & waterYear != "2021" & waterYear != "2022") & waterYear != "2014") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_839_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

ggplot(snotel_839_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_839_xts <- xts(snotel_839_clean_culled$temperature_mean, order.by = snotel_839_clean_culled$Date)

dygraph(temp_839_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
#snotel_839_clean_culled <- snotel_839_clean_culled %>% 
#  filter(temperature_mean < 19.3)

temp_839_xts <- xts(snotel_839_clean_culled$temperature_mean, order.by = snotel_839_clean_culled$Date)

dygraph(temp_839_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Upper Rio Grande 839

Oyler -> Morrisey 5/26/2004

snotel_839_adjusted <- snotel_839_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2004-05-26", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

839 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_839 <- snotel_839_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_839 <- yearly_wy_aver_839 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_839 <- daily_wy_aver_839 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_839$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_839 <-daily_wy_aver_839 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_839$date_temp <- signif(daily_wy_aver2_839$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_839, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

839 SD

standard_dev_839 <- daily_wy_aver_839 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_839 <- standard_dev_839 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_839 <- standard_dev_all_839 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_839 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 9.000013
1988 9.231875
1989 9.553651
1990 8.539684
1991 9.095321
1992 8.840589
1993 9.034837
1995 8.477482
1996 8.374844
1997 8.635328
1998 9.151866
1999 7.826807
2000 8.411485
2001 9.307564
2002 9.394376
2003 8.653297
2005 8.284320
2006 7.873798
2007 8.375674
2008 8.874800
2009 7.572710
2010 9.043721
2011 8.610385
2013 9.013859
2015 7.623354
2017 7.852672
2018 7.363181
2019 8.357738
2020 8.830386
2021 8.642675
2022 7.975856
ggplot(standard_dev_all_839, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 839 average temperatures for water years 2005-2021

MK & SS for 839 (non-corrected)

sd_mk_839 <- mk.test(standard_dev_all_839$sd_2)
print(sd_mk_839)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_839$sd_2
## z = -2.5495, n = 31, p-value = 0.01079
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -151.0000000 3461.6666667   -0.3247312
sd_sens_839 <- sens.slope(standard_dev_all_839$sd_2)
print(sd_sens_839)
## 
##  Sen's slope
## 
## data:  standard_dev_all_839$sd_2
## z = -2.5495, n = 31, p-value = 0.01079
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.051315272 -0.008246409
## sample estimates:
## Sen's slope 
## -0.03122783

Corrected

839 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_839_ad <- snotel_839_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_839_ad <- yearly_wy_aver_839_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_839_ad <- daily_wy_aver_839_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_839_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_839_ad <-daily_wy_aver_839_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_839_ad$date_temp_ad <- signif(daily_wy_aver2_839_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_839_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

839 SS (corrected)

standard_dev_839_ad <- daily_wy_aver_839_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_839_ad <- standard_dev_839_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_839_ad <- standard_dev_all_839_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_839_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 8.404486
1988 8.657820
1989 8.999163
1990 7.974877
1991 8.538551
1992 8.275130
1993 8.430766
1995 7.896968
1996 7.767745
1997 8.045426
1998 8.511594
1999 7.265049
2000 7.783278
2001 8.666438
2002 8.721243
2003 8.026530
2005 8.284037
2006 7.873931
2007 8.375662
2008 8.875360
2009 7.573403
2010 9.044246
2011 8.610915
2013 9.013498
2015 7.624141
2017 7.852753
2018 7.363589
2019 8.358178
2020 8.831049
2021 8.643116
2022 7.976653
ggplot(standard_dev_all_839_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 839 average temperatures for water years 1986-2021

MK & SS 839 (corrected)

sd_mk_839_ad <- mk.test(standard_dev_all_839_ad$sd_2)
print(sd_mk_839_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_839_ad$sd_2
## z = -0.16996, n = 31, p-value = 0.865
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##  -11.00000000 3461.66666667   -0.02365591
sd_sens_839_ad <- sens.slope(standard_dev_all_839_ad$sd_2)
print(sd_sens_839_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_839_ad$sd_2
## z = -0.16996, n = 31, p-value = 0.865
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.02519708  0.02245395
## sample estimates:
##  Sen's slope 
## -0.002233571

Upper San Juan 840

Morrisey 4/7/2004

snotel_840 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "840")
#str(snotel_840) # check the date, usually a character.  

snotel_840$Date <- as.Date(snotel_840$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_840_clean <- snotel_840 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_840_clean <- snotel_840_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_840_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_840_clean <- snotel_840_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40) %>% 
  filter(temperature_mean < 25)
ggplot(snotel_840_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
#snotel_840_cull_count <- snotel_840_clean %>% 
#  filter(temperature_min > -40) %>% 
#  count(waterYear)

#snotel_840_cull_count

# filtering for too few observations in a year
snotel_840_cull_count_days <- snotel_840_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_840_cull_count_days
## # A tibble: 11 x 2
## # Groups:   waterYear [11]
##    waterYear     n
##        <dbl> <int>
##  1      1980   310
##  2      1981   274
##  3      1982     1
##  4      1983   318
##  5      1984   322
##  6      1985   275
##  7      1993   349
##  8      1994   345
##  9      2005   309
## 10      2006   325
## 11      2013   314

1980, 1981, 1982, 1983, 1984, 1985, 1993, 1994, 2005, 2006, 2013 need to be culled.

snotel_840_clean_culled <- snotel_840_clean %>% 
  filter(waterYear >= "1986" & waterYear != "1993" & waterYear != "1994" & waterYear != "2005" & waterYear != "2006" & waterYear != "2013")# & waterYear != "2022") & waterYear != "2014") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_840_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

ggplot(snotel_840_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_840_xts <- xts(snotel_840_clean_culled$temperature_mean, order.by = snotel_840_clean_culled$Date)

dygraph(temp_840_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
#snotel_840_clean_culled <- snotel_840_clean_culled %>% 
#  filter(temperature_mean < 19.3)

temp_840_xts <- xts(snotel_840_clean_culled$temperature_mean, order.by = snotel_840_clean_culled$Date)

dygraph(temp_840_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Upper San Juan 840

Morrisey 4/7/2004

snotel_840_adjusted <- snotel_840_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2004-04-07", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

840 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_840 <- snotel_840_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_840 <- yearly_wy_aver_840 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_840 <- daily_wy_aver_840 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_840$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_840 <-daily_wy_aver_840 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_840$date_temp <- signif(daily_wy_aver2_840$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_840, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

840 SD

standard_dev_840 <- daily_wy_aver_840 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_840 <- standard_dev_840 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_840 <- standard_dev_all_840 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_840 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1986 7.232619
1987 7.739156
1988 8.527104
1989 8.409877
1990 7.857130
1991 8.153466
1992 7.489499
1995 7.461187
1996 7.863786
1997 7.870552
1998 8.397044
1999 7.262936
2000 7.859027
2001 8.336002
2002 8.619541
2003 8.032941
2004 8.221773
2007 7.731123
2008 7.978427
2009 7.095790
2010 8.100698
2011 7.956902
2012 7.805592
2014 7.587053
2015 6.894972
2016 7.701017
2017 7.249233
2018 7.180902
2019 7.751539
2020 7.888215
2021 7.886812
2022 7.530559
ggplot(standard_dev_all_840, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 840 average temperatures for water years 2005-2021

MK & SS for 840 (non-corrected)

sd_mk_840 <- mk.test(standard_dev_all_840$sd_2)
print(sd_mk_840)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_840$sd_2
## z = -1.5406, n = 32, p-value = 0.1234
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  -96.0000000 3802.6666667   -0.1935484
sd_sens_840 <- sens.slope(standard_dev_all_840$sd_2)
print(sd_sens_840)
## 
##  Sen's slope
## 
## data:  standard_dev_all_840$sd_2
## z = -1.5406, n = 32, p-value = 0.1234
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.031404215  0.001716915
## sample estimates:
## Sen's slope 
##  -0.0141416

Corrected

840 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_840_ad <- snotel_840_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_840_ad <- yearly_wy_aver_840_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_840_ad <- daily_wy_aver_840_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_840_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_840_ad <-daily_wy_aver_840_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_840_ad$date_temp_ad <- signif(daily_wy_aver2_840_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_840_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

840 SS (corrected)

standard_dev_840_ad <- daily_wy_aver_840_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_840_ad <- standard_dev_840_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_840_ad <- standard_dev_all_840_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_840_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1986 6.729334
1987 7.211134
1988 7.975188
1989 7.874323
1990 7.334504
1991 7.648553
1992 6.993487
1995 6.945265
1996 7.326356
1997 7.346542
1998 7.819710
1999 6.756828
2000 7.298996
2001 7.766886
2002 8.021992
2003 7.460019
2004 7.506468
2007 7.729865
2008 7.977294
2009 7.094470
2010 8.099982
2011 7.955897
2012 7.804828
2014 7.585948
2015 6.894260
2016 7.700243
2017 7.247788
2018 7.180489
2019 7.751133
2020 7.887037
2021 7.885810
2022 7.530352
ggplot(standard_dev_all_840_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 840 average temperatures for water years 1986-2021

MK & SS 840 (corrected)

sd_mk_840_ad <- mk.test(standard_dev_all_840_ad$sd_2)
print(sd_mk_840_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_840_ad$sd_2
## z = 1.2811, n = 32, p-value = 0.2002
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   80.0000000 3802.6666667    0.1612903
sd_sens_840_ad <- sens.slope(standard_dev_all_840_ad$sd_2)
print(sd_sens_840_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_840_ad$sd_2
## z = 1.2811, n = 32, p-value = 0.2002
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.005553565  0.027024754
## sample estimates:
## Sen's slope 
##   0.0120862

Vallecito 843 Morrisey 7/22/2005

snotel_843 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "843")
#str(snotel_843) # check the date, usually a character.  

snotel_843$Date <- as.Date(snotel_843$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_843_clean <- snotel_843 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_843_clean <- snotel_843_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_843_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_843_clean <- snotel_843_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40) %>% 
  filter(temperature_mean < 25)
ggplot(snotel_843_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
#snotel_843_cull_count <- snotel_843_clean %>% 
#  filter(temperature_min > -40) %>% 
#  count(waterYear)

#snotel_843_cull_count

# filtering for too few observations in a year
snotel_843_cull_count_days <- snotel_843_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_843_cull_count_days
## # A tibble: 2 x 2
## # Groups:   waterYear [2]
##   waterYear     n
##       <dbl> <int>
## 1      1994   339
## 2      2003   313

1994 & 2003 need to be culled.

snotel_843_clean_culled <- snotel_843_clean %>% 
  filter(waterYear != "1994" & waterYear != "2003")# & waterYear != "1994" & waterYear != "2005" & waterYear != "2006" & waterYear != "2013")# & waterYear != "2022") & waterYear != "2014") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_843_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

ggplot(snotel_843_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_843_xts <- xts(snotel_843_clean_culled$temperature_mean, order.by = snotel_843_clean_culled$Date)

dygraph(temp_843_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
#snotel_843_clean_culled <- snotel_843_clean_culled %>% 
#  filter(temperature_mean < 19.3)

temp_843_xts <- xts(snotel_843_clean_culled$temperature_mean, order.by = snotel_843_clean_culled$Date)

dygraph(temp_843_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Vallecito 843

Morrisey 7/22/2005

snotel_843_adjusted <- snotel_843_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2005-07-22", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

843 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_843 <- snotel_843_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_843 <- yearly_wy_aver_843 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_843 <- daily_wy_aver_843 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_843$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_843 <-daily_wy_aver_843 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_843$date_temp <- signif(daily_wy_aver2_843$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_843, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

843 SD

standard_dev_843 <- daily_wy_aver_843 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_843 <- standard_dev_843 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_843 <- standard_dev_all_843 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_843 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 7.744307
1988 8.133975
1989 8.353802
1990 7.814664
1991 7.913458
1992 7.258470
1993 8.017543
1995 7.737939
1996 8.012757
1997 8.122829
1998 8.284299
1999 6.874497
2000 8.036954
2001 8.326535
2002 8.475624
2004 5.169906
2005 7.879202
2006 7.087987
2007 7.572592
2008 8.035842
2009 7.216959
2010 7.825626
2011 7.886538
2012 7.648815
2013 8.040098
2014 7.292347
2015 6.734206
2016 7.523001
2017 7.419976
2018 7.106676
2019 7.826098
2020 7.899858
2021 7.970432
2022 7.307690
ggplot(standard_dev_all_843, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 843 average temperatures for water years 2005-2021

MK & SS for 843 (non-corrected)

sd_mk_843 <- mk.test(standard_dev_all_843$sd_2)
print(sd_mk_843)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_843$sd_2
## z = -1.6307, n = 34, p-value = 0.103
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##           S        varS         tau 
## -111.000000 4550.333333   -0.197861
sd_sens_843 <- sens.slope(standard_dev_all_843$sd_2)
print(sd_sens_843)
## 
##  Sen's slope
## 
## data:  standard_dev_all_843$sd_2
## z = -1.6307, n = 34, p-value = 0.103
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.029639051  0.001693822
## sample estimates:
## Sen's slope 
## -0.01383409

Corrected

843 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_843_ad <- snotel_843_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_843_ad <- yearly_wy_aver_843_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_843_ad <- daily_wy_aver_843_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_843_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_843_ad <-daily_wy_aver_843_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_843_ad$date_temp_ad <- signif(daily_wy_aver2_843_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_843_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

843 SS (corrected)

standard_dev_843_ad <- daily_wy_aver_843_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_843_ad <- standard_dev_843_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_843_ad <- standard_dev_all_843_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_843_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 7.167709
1988 7.554905
1989 7.775142
1990 7.236968
1991 7.363376
1992 6.713681
1993 7.440691
1995 7.144630
1996 7.409097
1997 7.529924
1998 7.661187
1999 6.355448
2000 7.405969
2001 7.695693
2002 7.828232
2004 4.700776
2005 7.205286
2006 7.087285
2007 7.571307
2008 8.034938
2009 7.215649
2010 7.824301
2011 7.885342
2012 7.647841
2013 8.038906
2014 7.290741
2015 6.733530
2016 7.521405
2017 7.418681
2018 7.106001
2019 7.824924
2020 7.898803
2021 7.969378
2022 7.306728
ggplot(standard_dev_all_843_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 843 average temperatures for water years 1986-2021

MK & SS 843 (corrected)

sd_mk_843_ad <- mk.test(standard_dev_all_843_ad$sd_2)
print(sd_mk_843_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_843_ad$sd_2
## z = 1.601, n = 34, p-value = 0.1094
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  109.0000000 4550.3333333    0.1942959
sd_sens_843_ad <- sens.slope(standard_dev_all_843_ad$sd_2)
print(sd_sens_843_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_843_ad$sd_2
## z = 1.601, n = 34, p-value = 0.1094
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.003122605  0.025938612
## sample estimates:
## Sen's slope 
##  0.01173229

Wolf Creek Summit 874

Morrisey 7/12/2004

snotel_874 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "874")
#str(snotel_874) # check the date, usually a character.  

snotel_874$Date <- as.Date(snotel_874$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_874_clean <- snotel_874 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_874_clean <- snotel_874_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_874_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_874_clean <- snotel_874_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40) %>% 
  filter(temperature_mean < 25)
ggplot(snotel_874_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
#snotel_874_cull_count <- snotel_874_clean %>% 
#  filter(temperature_min > -40) %>% 
#  count(waterYear)

#snotel_874_cull_count

# filtering for too few observations in a year
snotel_874_cull_count_days <- snotel_874_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_874_cull_count_days
## # A tibble: 11 x 2
## # Groups:   waterYear [11]
##    waterYear     n
##        <dbl> <int>
##  1      1990   295
##  2      1994   313
##  3      2003   311
##  4      2006   340
##  5      2007   340
##  6      2008   315
##  7      2009   320
##  8      2010   303
##  9      2013   304
## 10      2016   338
## 11      2017   314

1990, 1994, 2003, 2006, 2007, 2008, 2009, 2010, 2013, 2016, 2017 need to be culled.

snotel_874_clean_culled <- snotel_874_clean %>% 
  filter(waterYear != "1990" & waterYear != "1994" & waterYear != "2003" & waterYear != "2006" & waterYear != "2007" & waterYear != "2008" & waterYear != "2009" & waterYear != "2010" & waterYear != "2013" & waterYear != "2016" & waterYear != "2017") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_874_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

ggplot(snotel_874_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_874_xts <- xts(snotel_874_clean_culled$temperature_mean, order.by = snotel_874_clean_culled$Date)

dygraph(temp_874_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
#snotel_874_clean_culled <- snotel_874_clean_culled %>% 
#  filter(temperature_mean < 19.3)

temp_874_xts <- xts(snotel_874_clean_culled$temperature_mean, order.by = snotel_874_clean_culled$Date)

dygraph(temp_874_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Wolf Creek Summit 874

Morrisey 7/12/2004

snotel_874_adjusted <- snotel_874_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2004-07-12", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

874 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_874 <- snotel_874_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_874 <- yearly_wy_aver_874 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_874 <- daily_wy_aver_874 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_874$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_874 <-daily_wy_aver_874 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_874$date_temp <- signif(daily_wy_aver2_874$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_874, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

874 SD

standard_dev_874 <- daily_wy_aver_874 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_874 <- standard_dev_874 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_874 <- standard_dev_all_874 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_874 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 7.854854
1988 8.380441
1989 8.376871
1991 7.952272
1992 7.274045
1993 7.924599
1995 7.536845
1996 8.020204
1997 8.005351
1998 8.406529
1999 7.057906
2000 8.042917
2001 8.360075
2002 8.824249
2004 8.065839
2005 7.482529
2011 8.051681
2012 7.789860
2014 7.447867
2015 6.928269
2018 7.223986
2019 7.800112
2020 7.863897
2021 8.007803
2022 7.236522
ggplot(standard_dev_all_874, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 874 average temperatures for water years 2005-2021

MK & SS for 874 (non-corrected)

sd_mk_874 <- mk.test(standard_dev_all_874$sd_2)
print(sd_mk_874)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_874$sd_2
## z = -1.4247, n = 25, p-value = 0.1543
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  -62.0000000 1833.3333333   -0.2066667
sd_sens_874 <- sens.slope(standard_dev_all_874$sd_2)
print(sd_sens_874)
## 
##  Sen's slope
## 
## data:  standard_dev_all_874$sd_2
## z = -1.4247, n = 25, p-value = 0.1543
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.046813576  0.005791174
## sample estimates:
## Sen's slope 
## -0.01681136

Corrected

874 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_874_ad <- snotel_874_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_874_ad <- yearly_wy_aver_874_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_874_ad <- daily_wy_aver_874_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_874_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_874_ad <-daily_wy_aver_874_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_874_ad$date_temp_ad <- signif(daily_wy_aver2_874_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_874_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

874 SS (corrected)

standard_dev_874_ad <- daily_wy_aver_874_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_874_ad <- standard_dev_874_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_874_ad <- standard_dev_all_874_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_874_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 7.327962
1988 7.851838
1989 7.848099
1991 7.458480
1992 6.775174
1993 7.395057
1995 7.010946
1996 7.470539
1997 7.481361
1998 7.826754
1999 6.568878
2000 7.465159
2001 7.779500
2002 8.206390
2004 7.444487
2005 7.482090
2011 8.052315
2012 7.790524
2014 7.447936
2015 6.928554
2018 7.224892
2019 7.800202
2020 7.864350
2021 8.008044
2022 7.236992
ggplot(standard_dev_all_874_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 874 average temperatures for water years 1986-2021

MK & SS 874 (corrected)

sd_mk_874_ad <- mk.test(standard_dev_all_874_ad$sd_2)
print(sd_mk_874_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_874_ad$sd_2
## z = 0.81742, n = 25, p-value = 0.4137
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##        S     varS      tau 
##   36.000 1833.333    0.120
sd_sens_874_ad <- sens.slope(standard_dev_all_874_ad$sd_2)
print(sd_sens_874_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_874_ad$sd_2
## z = 0.81742, n = 25, p-value = 0.4137
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.006890724  0.033504267
## sample estimates:
## Sen's slope 
## 0.006882736

Figure captions are not correct.